“ R is a free software environment for statistical computing and graphics. It compiles and runs on a wide variety of UNIX platforms, Windows and MacOS. ”
R의 기본 연산 단위는 벡터이다.
x=c(1,2,3,4,5,6) ## vector of variable
y=c(7,8,9,10,11,12)
x+y; x*y; sqrt(x)
## [1] 8 10 12 14 16 18
## [1] 7 16 27 40 55 72
## [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490
mean(x); sd(x); length(x)
## [1] 3.5
## [1] 1.870829
## [1] 6
x[2]; x[-2]; x[c(1,3,5)]
## [1] 2
## [1] 1 3 4 5 6
## [1] 1 3 5
벡터를 생성하는 다양한 방법
## Sequence
v1=seq(-5,5,by=.2); v1
## [1] -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2
## [16] -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8
## [31] 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8
## [46] 4.0 4.2 4.4 4.6 4.8 5.0
## Repeat
v2=rep(1,3); v2
## [1] 1 1 1
v3=rep(c(1,2,3),2); v3 ## Repeat for vector
## [1] 1 2 3 1 2 3
v4=rep(c(1,2,3),each = 2); v4 ## Repeat for vector : each
## [1] 1 1 2 2 3 3
## for loop
for (i in 1:3){
print(i)
}
## [1] 1
## [1] 2
## [1] 3
i=0
for (j in c(1,2,4,5,6)){
i=i+j
}
i
## [1] 18
## if
x=5
if (x >=3 ){
x=x+3
}
x
## [1] 8
x=5
if (x >=10){
print("High")
} else if (x >=5){
print("Medium")
} else {
print("Low")
} ## if, else if 주의: 반드시 } 와 같은 줄에 위치하도록.
## [1] "Medium"
## ifelse
x=5
y=ifelse(x==5,"OK","Suck") ## ifelse(조건,참일때,거짓일때)
y
## [1] "OK"
R은 벡터 기반의 연산을 지원하므로 이를 이용하면 여러가지 계산을 간단히 할 수 있다.
apply는 행 또는 열 단위의 연산을 쉽게 할 수 있도록 지원하는 함수이다.
mat=matrix(1:20,nrow=4,byrow=T) ## 4행 5열, byrow=T : 행부터 채운다.
mat
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 2 3 4 5
## [2,] 6 7 8 9 10
## [3,] 11 12 13 14 15
## [4,] 16 17 18 19 20
apply(mat,1,mean) ## 1: 행- 행별로 평균
## [1] 3 8 13 18
apply(mat,2,mean) ## 2: 열- 열별로 평균
## [1] 8.5 9.5 10.5 11.5 12.5
sapply는 실행결과를 벡터로 출력하고 lapply는 실행결과를 리스트로 출력한다.
sapply(1:nrow(mat),function(x){mean(mat[x,])})
## [1] 3 8 13 18
sapply(1:ncol(mat),function(x){mean(mat[,x])})
## [1] 8.5 9.5 10.5 11.5 12.5
lapply(1:nrow(mat),function(x){mean(mat[x,])})
## [[1]]
## [1] 3
##
## [[2]]
## [1] 8
##
## [[3]]
## [1] 13
##
## [[4]]
## [1] 18
unlist(lapply(1:nrow(mat),function(x){mean(mat[x,])})) ## Same to sapply
## [1] 3 8 13 18
for문의 과도한 사용은 분석의 속도를 저하시키는 가장 큰 요인 중 하나이다.
따라서 for문을 최대한 덜 쓰고 벡터연산을 활용하는 것이 실행시간을 줄이는 핵심이다.
## for 문을 이용한 합구하는 함수 만들기
sum_f=function(x){
out=0
for (i in 1:x){
out=out+i
}
return(out)
}
system.time(sum_f(10^6)) ## system.time - 시간 잰다.
## user system elapsed
## 0.048 0.000 0.048
system.time(sum(as.numeric(1:10^6))) ## 내장 sum 함수- 더할 것이 많으면 as.numeric을 넣어줘야 한다.
## user system elapsed
## 0 0 0
데이터를 불러오기 전에 미리 디렉토리를 지정하면 그 다음부터는 편하게 읽고 쓸 수 있다.
여기서 주의해야 할 점은 폴더간의 구분을 / 로 해야 한다는 점이다. R 은 유닉스 기반의 언어이기 때문이다.
setwd("/home/heejooko/ShinyApps/introduction to R")
getwd()
## [1] "/home/heejooko/ShinyApps/introduction to R"
예제로 파일을 불러와 살펴보자.
a<-read.csv("practice.csv")
a
head(a); tail(a)
a$SEX <- ifelse(a$SEX=='0',NA,a$SEX)
tail(a)
class(a)
## [1] "data.frame"
names(a)
## [1] "CASEID" "AGE" "SEX" "HEIGHT" "WEIGHT" "DM"
## [7] "HTN" "PH" "TREATMENT"
nrow(a); ncol(a)
## [1] 50
## [1] 9
a[,1:5]
a[20,]
magrittr 패키지의 %>% 연산자를 이용하면 직관적인 코드를 작성할 수 있다.
# install.packages("magrittr")
library(magrittr)
a %>% head ## head(a)와 같다.
a %>% head %>% nrow ## nrow(head(a))와 같다.
## [1] 6
키에 대한 통계량을 살펴보자
mean(a$HEIGHT) ## 결측치 없을 때만
## [1] NA
mean(a$HEIGHT, na.rm=T) ## 결측치 빼고
## [1] 162.6043
round(mean(a$HEIGHT, na.rm=T), 1) ## 반올림
## [1] 162.6
a$HEIGHT %>% mean(.,na.rm=T) %>% round(.,1) ## 위와 같은 명령
## [1] 162.6
아래의 함수 모두 na.rm=T을 빠트려서는 안된다.
sd(a$HEIGHT,na.rm=T) ## 표준편차
## [1] 18.66815
var(a$HEIGHT,na.rm=T) ## 분산
## [1] 348.5
median(a$HEIGHT,na.rm=T) ## 중간값
## [1] 163.9
IQR(a$HEIGHT,na.rm=T) ## 25%-75% range
## [1] 14.1
quantile(a$HEIGHT, na.rm=T) ## quantile
## 0% 25% 50% 75% 100%
## 51.3 158.4 163.9 172.5 183.0
max(a$HEIGHT,na.rm=T) ## Max
## [1] 183
min(a$HEIGHT,na.rm=T) ## Min
## [1] 51.3
HEIGHT, WEIGHT 변수를 이용해서 BMI 변수를 만들어보자.
a$BMI <- a$WEIGHT/((a$HEIGHT/100)^2)
a$BMI
## [1] 29.18599 23.97353 NA 19.08418 22.53995 27.05627 22.75956
## [8] 22.72571 23.40672 22.32128 NA 25.27620 22.33710 19.48452
## [15] 23.63281 45.48324 34.88653 27.16967 18.36706 16.09042 22.35622
## [22] 20.16475 22.31107 21.05703 189.23201 21.54913 22.65067 16.23245
## [29] 22.75374 35.71304 NA 12.69589 24.13891 19.75689 20.04338
## [36] 19.21222 20.67644 29.01123 23.17902 NA 21.73339 16.15029
## [43] 21.01780 23.50206 NA 23.80408 20.75674 24.75186 28.06813
## [50] 17.61528
25를 기준으로 BMI 분류 변수를 만들어보자
a$BMI_cat=0 ## 0으로 된 변수만들고
a$BMI_cat[a$BMI>=25]=1 ## 25이상인 것만 1로 바꾼다.
table(a$BMI_cat)
##
## 0 1
## 40 10
a$BMI_cat=a$BMI>=25 ## 방법 2: TRUE of FALSE
table(a$BMI_cat)
##
## FALSE TRUE
## 35 10
a$BMI_cat=ifelse(a$BMI>=25,"1","0") ## 방법 3: ifelse문 이용
table(a$BMI_cat)
##
## 0 1
## 35 10
ifelse의 응용
a$BMI_cat=ifelse(a$BMI>=25,"overweight",ifelse(a$BMI<18.5,"underweight","normal"))
table(a$BMI_cat)
##
## normal overweight underweight
## 29 10 6
a
a 자료의 대략적인 요약을 살펴보자.
summary(a)
## CASEID AGE SEX HEIGHT
## Min. : 1.00 Min. :15.00 Length:50 Min. : 51.3
## 1st Qu.:13.25 1st Qu.:27.50 Class :character 1st Qu.:158.4
## Median :25.50 Median :46.50 Mode :character Median :163.9
## Mean :25.50 Mean :44.32 Mean :162.6
## 3rd Qu.:37.75 3rd Qu.:59.50 3rd Qu.:172.5
## Max. :50.00 Max. :75.00 Max. :183.0
## NA's :3
## WEIGHT DM HTN PH TREATMENT
## Min. : 40.00 Min. :0.00 Min. :0.00 Min. :6.850 Min. :0.0
## 1st Qu.: 54.75 1st Qu.:0.00 1st Qu.:0.00 1st Qu.:7.122 1st Qu.:0.0
## Median : 60.00 Median :1.00 Median :0.00 Median :7.180 Median :0.5
## Mean : 62.78 Mean :0.66 Mean :0.46 Mean :7.176 Mean :0.5
## 3rd Qu.: 67.58 3rd Qu.:1.00 3rd Qu.:1.00 3rd Qu.:7.240 3rd Qu.:1.0
## Max. :120.40 Max. :1.00 Max. :1.00 Max. :7.370 Max. :1.0
## NA's :4
## BMI BMI_cat
## Min. : 12.70 Length:50
## 1st Qu.: 20.16 Class :character
## Median : 22.54 Mode :character
## Mean : 26.80
## 3rd Qu.: 24.14
## Max. :189.23
## NA's :5
이상한 점은 변수의 class가 올바르게 설정되지 않았기 때문이다.
sapply(a,class)
## CASEID AGE SEX HEIGHT WEIGHT DM
## "integer" "integer" "character" "numeric" "numeric" "integer"
## HTN PH TREATMENT BMI BMI_cat
## "integer" "numeric" "integer" "numeric" "character"
CASEID는 character로, SEX, DM, HTN, TREATMENT, BMI_cat 변수는 범주형 변수로 바꾸고자 한다.
a$CASEID <- as.character(a$CASEID)
cat_vars <- c("SEX","DM","HTN","TREATMENT","BMI_cat")
for(vn in cat_vars){
a[,vn] <- as.factor(a[,vn])
}
sapply(a,class)
## CASEID AGE SEX HEIGHT WEIGHT DM
## "character" "integer" "factor" "numeric" "numeric" "factor"
## HTN PH TREATMENT BMI BMI_cat
## "factor" "numeric" "factor" "numeric" "factor"
summary(a)
## CASEID AGE SEX HEIGHT WEIGHT
## Length:50 Min. :15.00 F :19 Min. : 51.3 Min. : 40.00
## Class :character 1st Qu.:27.50 M :29 1st Qu.:158.4 1st Qu.: 54.75
## Mode :character Median :46.50 NA's: 2 Median :163.9 Median : 60.00
## Mean :44.32 Mean :162.6 Mean : 62.78
## 3rd Qu.:59.50 3rd Qu.:172.5 3rd Qu.: 67.58
## Max. :75.00 Max. :183.0 Max. :120.40
## NA's :3 NA's :4
## DM HTN PH TREATMENT BMI BMI_cat
## 0:17 0:27 Min. :6.850 0:25 Min. : 12.70 normal :29
## 1:33 1:23 1st Qu.:7.122 1:25 1st Qu.: 20.16 overweight :10
## Median :7.180 Median : 22.54 underweight: 6
## Mean :7.176 Mean : 26.80 NA's : 5
## 3rd Qu.:7.240 3rd Qu.: 24.14
## Max. :7.370 Max. :189.23
## NA's :5
특정 조건을 만족하는 부분 데이터셋을 만드는 것은 지금까지 배웠던 것을 응용할 수도 있고 subset이라는 함수를 쓸 수도 있다.
a[a$AGE>=60 & a$SEX=='F',] ## 행에 조건 써준다.
subset(a,a$AGE>=60 & a$SEX=='F') ## 문법이 간단하여 이것을 주로 쓴다.
subset(a,a$AGE>=60 & a$SEX=='F', select=c("AGE","SEX","DM")) ## select로 원하는 변수만 추려낸다.
대부분의 의학 논문에는 데이터의 특성을 나타내는 Table1이 포함되어야 한다.
Table1을 만드는 다양한 함수 중 tableone 패키지의 CreateTableOne 함수를 사용해보겠다.
# install.packages("tableone")
library(tableone)
tb1=CreateTableOne(vars = names(a[,2:9]), ## 테이블에 들어갈 변수들
strata = "BMI_cat", ## 그룹 변수(여러 개 가능)
factorVars = cat_vars, ## 범주형 변수들
data = a,
includeNA = T ## 범주형 변수에서 NA를 하나의 범주로 포함할 것인가?
)
tb1
## Stratified by BMI_cat
## normal overweight underweight p test
## n 29 10 6
## AGE (mean (SD)) 46.17 (15.78) 43.10 (22.98) 42.00 (21.12) 0.823
## SEX (%) 0.805
## F 9 (31.0) 5 (50.0) 2 (33.3)
## M 19 (65.5) 5 (50.0) 4 (66.7)
## NA 1 ( 3.4) 0 ( 0.0) 0 ( 0.0)
## HEIGHT (mean (SD)) 165.07 (8.25) 148.58 (34.76) 175.68 (2.45) 0.010
## WEIGHT (mean (SD)) 59.98 (6.88) 76.81 (20.08) 49.96 (6.00) <0.001
## DM = 1 (%) 20 (69.0) 5 (50.0) 4 (66.7) 0.554
## HTN = 1 (%) 17 (58.6) 3 (30.0) 2 (33.3) 0.211
## PH (mean (SD)) 7.17 (0.10) 7.21 (0.07) 7.17 (0.18) 0.651
## TREATMENT = 1 (%) 15 (51.7) 4 (40.0) 2 (33.3) 0.636
R의 큰 장점 중 하나는 수려한 그래픽이다. 간단히 예제 데이터로 두 가지를 그려보자.
오늘은 데이터의 형태에 따라 상상하는 거의 모든 그래프를 그릴 수 있다는 것 정도만 기억하자.
library(ggplot2)
a<-na.omit(a) ## 결측치가 없는 행만 추린다
ggplot(a, aes(x=HEIGHT, y=BMI, color=SEX)) +
geom_point() +
geom_smooth(method=lm, se=FALSE, fullrange=TRUE)
## `geom_smooth()` using formula 'y ~ x'
a$BMI_cat <- factor(a$BMI_cat, levels = c("underweight", "normal", "overweight"))
ggplot(a, aes(x=BMI_cat, y=PH, color=BMI_cat)) +
geom_boxplot() +
geom_jitter(color="black", size=0.4, alpha=0.9)
다양한 패키지, 함수, 옵션들을 활용하여 이런 표와 그림들을 그릴 수 있다.
그림을 잘 그리기 위해서는 데이터의 적절한 가공, 적절한 함수 선택, 가시성을 높여주는 적절한 옵션의 추가가 중요하다.
datatable 함수에 formatStyle로 디자인했다. Figure 3. Differences in dose distribution values
The values represent the differences between the values obtained by the techniques listed in the columns and those obtained by the techniques listed in the rows. Red indicates statistically significant-positive differences and blue indicates statistically significant-negative differences
ggplot, geom_line 함수를 사용했다.